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Abstract. We present a new algorithm based on a Cartesian mesh for the numerical ap- 
proximation of kinetic models for chemosensitive movements set in an arbitrary geometry. 
We investigate the influence of the geometry on the collective behavior of bacteria described 
by a kinetic equation interacting with nutrients and chemoattractants. Numerical simula- 
tions are performed to verify accuracy and stability of the scheme and its ability to exhibit 
aggregation of cells and wave propagations. Finally some comparisons with experiments 
show the robustness and accuracy of such kinetic models. 
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1. Introduction 

This paper is devoted to the numerical simulation of the process of cellular spatial orga- 
nization driven by chemotaxis. The effective mechanism by which individual cells undergo 
directed motion varies among organisms. Here we are particularly interested in bacterial 
migration, characterized by the smallness of the cells, and their ability to move to several 
orders of magnitude in the attractant concentration. Several models, depending on the level 
of description, have been developed mathematically for the collective motion of cells. We 
refer to [5, 24, 25] for a review on parabolic, hyperbolic and kinetic models and to [^^ '^^^ 28] 
for traveling waves drivn by growth and chemotaxis. Among them the kinetic model in- 
troduced by Othmer, Dunbar and Alt [ , ], describes a population of bacteria in motion 
(for instance the E. Coli) in interactions with a chemoattractant concentration [ ]. These 
cells are so small that they are not able to measure any head-to-tail gradient of the chemical 
concentration, and to choose directly some preferred direction of motion towards high con- 
centrated regions. Therefore they develop an indirect strategy to select favorable regions, by 
detecting a sort of time derivative in the concentration along their pathways, and to react 
accordingly [16]. 

More precisely, they undergo a jump process where free runs are followed by a reorientation 
phenomenon called tumble [32]. For instance it is known that E. Coli increases the time 
spent in running along a favourable direction [ ]. This jump process can be described by two 
different informations. First cells switch the rotation of their flagella, from counter-clockwise. 

The authors are partially supported by the European Research Council ERC Starting Grant 2009, project 
2?>9m?>-NuSiKiMo. 
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called free runs, to clockwise called reorientation or tumbling phase, and conversely. This 
decision is the result of a complex chain of reactions inside the cells, driven by the external 
concentration of the chemoattractant [ ^9, 32]. Then bacteria select a new direction, but 
they are unable to choose directly a favourable direction, so they randomly choose a new 
orientation. During the "run" phases a bacterium moves with a constant speed in a given 
direction while during a "tumbling" event it changes direction randomly. 

2. Kinetic models for bacterial chemotaxis 

In the simple situation, C. S. Patlak [23] and E. F. Keller & L. A. Segel [20] considered 
a density of cells which interacts with two chemical substances. The cells consume nutrients 
which drive the migration and excrete a chemoattractant that prevents the dispersion of 
the population. However, this approach is not always sufficiently precise to describe the 
evolution of bacteria movements. Hence, this phenomenon of run and tumble can be modeled 
by a stochastic process called the velocity-jump process, and has been introduced by Alt [^] 
and further developed in [ ]. A kinetic transport model to describe this velocity jump 
process consists to study the evolution of the bacterial population by the local density of cells 
/(t,x, v) > located in position x, at time t and swimming in the direction v g y. 

Here the set V of all possible velocities is bounded and symmetric in general. In two 
dimensions, the modulus of the speed is a constant vq, hence V = S{0^vo) circle centred in 
with a radius vq > 0. 

A kinetic transport model to describe this velocity-jump process has been introduced by 
W. Alt [ ] inspired by the Boltzmann equation [22] where the tumbles appear as scattering 
events and all the fluxes are explicitly introduced [ ]. We consider the Boltzmann type 
equation 

(2.1) ^ + vVx/ = g(/)+r/, x€0, ^r€V, 
where Q{f) is the Boltzmann type tumbling operator 

Qif) = r(v, v') fit, X, v') dv' - Tiv', v) fit, X, v) dv'. 

The contribution of the tumbles is introduced with a transition (scattering) kernel r(v, v^) > 
which stands for the change of velocity from to v; r is the division rate of the bacteria 
(r = In2/T2 where T2 is the mean doubling time). This equation is indeed a variant of the 
Boltzmann equation for gases, where collisions are delocalized via the secretion or consump- 
tion of chemical cues. 

In (2.1), the transition kernel T also depends on the local concentration of chemoattractant 
S'(t,x) and nutrient A^(t,x). To estimate the respective contributions on pulse speed of the 
biais of the run lengths and of preferential reorientation, it is possible to split this transition 
kernel T(v,v^) in two contributions, one being the tumbling rate A(v^), and the other one 
the reorientation effect during tumbles A^(v,v^): 

(2.2) r(v,v') = A(v')i^(v,v'), 

with the condition 

(2.3) J^K{wy)dw = 1, 

where the function K accounts for the persistence of the trajectories. For simplicity we 
consider the case of the absence of such an angular persistence, hence the turning kernel 
T(v,v^) is only proportional to the tumbling rate A(v^), i.e. K is constant. 

For the tumbling rate A(v^), we assume that bacteria are sensitive to the temporal vari- 
ations of attractant molecules via a logarithmic sensing mechanism [4, 15]. Therefore, the 
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tumble frequency only depends on the local gradients of nutrient and attract ant, both gradi- 
ents having independent and additive contributions. This gives 

A(v') = ^ (A^(v') + As(v')) 



(^^(atlog(iV) + V.Vxlog(7V)) ^ ^Ps {dtlog{S) + V.Vxlog(5))), 



2 
1 
2 

The nutrient and chemoattractant response functions i/jn and tps ^ire both positive and 
decreasing, expressing that cells are less likely to tumble (thus perform longer runs) when the 
external chemical signal increases. These functions are smooth and characterized by their 
characteristic time 6^ and 6^^ and their tumble frequency xa^ and xs- 

Here, we have chosen the following analytical form that encompassed these characteristics: 

(2.4) Va(^) = V'o(l-Xatanh(^)), a = {N,S}, 

where '0o is the mean tumbling frequency, and the parameters xn and xs are the modulation 
of tumble frequency. 

In order to define completely the mathematical problem (2.1), suitable boundary conditions 
on dQ should be applied. Here we consider wall type boundary conditions, for which emerging 
particles have been reflected elastically at the wall. More precisely, for x g d^} the smooth 
boundary dfl is assumed to have a unit inward normal n(x) and for v • n(x) > 0, we assume 
that at the solid boundary we have 

(2.5) /(t,x,v) = 7e[/(t,x,v)], ^€dn, vn(x)>0, 
with 

(2.6) 7^[/(t, X, v)] = fit, X, V - 2(v • n(x))n(x)). 

This boundary condition (2.5) guarantees the global conservation of mass [^]. 

The equations describing the behaviors of nutrients density N and chemottractant S are 
the same as in [ 1: 



dS 

(2.7) 



dt 
ON 



DsAS = -aS + b J^f{t,yi,v)dv, x g f]. 



DnAN = -cN / /(t,x,v)rfv, xGf7, 
^ ot Jv 

where a, b and c are respectively the degradation rate of the chemoattractant, its production 
rate and the consumption rate of the nutrient by the bacteria, whereas Ds and are the 
molecular diffusion coefficients. Finally, these equations are completed with homogeneous 
Neumann boundary conditions, i.e. 

(2.8) Vx«-n(x) = 0, a = {N,S}, yiedQ. 

The purpose of this work is to present a numerical scheme for (2.1)-(2.8) and to investigate 
numerically the occurrence of cells aggregation, pattern formation or travelling waves when 
it takes place, and the convergence to equilibrium otherwise for different geometries. Several 
numerical methods have already been developed to solve the Patlak-Keller-Segel model for 
chemotaxis using finite element methods [ ], finite volume methods [8, 9, 13], and the refer- 
ences therein. Other models have also been investigated numerically [11, 12, 21]. However, it 
seems that none of the above-mentioned numerical approaches have been studied for kinetic 
models (2.1)-(2.8). In the present paper we propose a numerical scheme for (2.1)-(2.8) and 
investigate the influence of the geometry on the collective behavior of bacteria. 
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We now briefly outline the contents of the paper. In the next section, we introduce the 
numerical approximation of (2.1) and (2.7) and describe the numerical approximation of 
the boundary (2.5), (2.6), (2.8). Two points are worth mentioning here. First, we restrict 
ourselves to the case of specular reflection which seems to be the most appropriate for the 
study of bacteria. One difficulty in the approximation of kinetic models for chemotaxis, is 
related to the fact that it can exhibits very different phenomena as finite time blow-up, cell 
aggregation, wave propagations. At the discrete level, our approximation should also be able 
to describe a similar property. Secondly, an important step is to discretize appropriately the 
effect of boundary. 

The final section is devoted to numerical simulations performed with the numerical scheme 
presented in Section 3. We investigate numerically cells aggregation, convergence to equilib- 
rium, and wave propagation in a bounded domain. 



3. Numerical resolution 
3.1. Numerical resolution of the kinetic model (2.1). We consider a computational 

domain [Xmin,Xniax] ^ [l/min? ymax] ^ ^5 SUCh that Q C [Xmin, Xmax] ^ [l/min? I/max] • 

The computational domain is covered by an uniform Cartesian mesh X/^ x V/^, where X/^, 
V/^ are defined by 



(3.1) 



. V/, := {v^- = vo{cos0j,smej), Oj = {j ^ 1/2) Av,0 < j < - 1} . 



The mesh steps are respectively Ax = (xmax-^min)/^x5 = (l/max-ymin)/^?/ and Av = 27v/ny. 

Let us denote fl^j an approximation of the distribution function /(t^, x^, Vj). We introduce 
the following finite difference scheme 

(3.2) ' ' ^ ^ ■ = ^'^^^^■^ ^ ■ 

where h - {At, Ax, Ay), v- V:sL,hfi'j ^ second-order approximation [ '] of the transport 
operator v- Vx/, and Qhif^j) approximation of the Boltzmann type tumbling operator 
Qif)' We will now focus on searching the approximation 
By using the trapezoidal rule, we have 



where {Th)^^-\^^ is an approximation of the transition kernel r(v,v^). We assume that the 
nutrient density N and the chemottractant S are known at time Moreover with the 

hypothesis that K is constant, the condition (2.3) implies that K = l/27r. Thus it remains to 
search an approximation of the tumbling rate, i.e. A^^^^^. It is also equivalent to search the 

local gradients of nutrient {Xn)^^^^^'^ and attractant (A^')^^''^^^. 

We study only the local gradients of nutrient (Aat)^^^^^, since the local gradient of at- 
tractant {Xs)^^^^^'^ has the same expression as (Aat)^^''^^^. We discretize the local gradient of 
nutrient as follows 

xn+l/2 1 , ^ ^jn+ll2\ 
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where D^N^^^^^ is a discrete time derivative and will be given in the section 3.2. Moreover 
we use centred difference approximation for • Vx^, which yields 

^i'V^,hN. ' =vocos9i ^+7;osin^^ ~^Ay 

In summary, the discrete tumbling rate reads 



.n+l/2 1//. \n+l/2] 



2 

Finally we reduce the tumbling operator as follows 

iO f -pn \ _ ^sr<^ ( A n+1/2 nn \ a n+1/2 nn 

3.2. Calculate the chemoattractant S or the nutrient N. The equations (2.7) for 
nutrient density N and chemottractant S are parabolic equations with source terms depending 
on the distribution function of density /. We study again the discretization for nutrient 
density, since the one for chemottractant is similar. The Euler implicit scheme is used for 
time integration. Hence the scheme for nutrient density N reads 

(3.3) D.N:"^''' - DM^hK^'l^ = -cK^'^'pl, 
where D^N^*^^'^ is an approximation of time derivative as follows 

(3.4) D,Nr'" = ^ ^ . 

Then we use a five points finite difference scheme to discretize /\N 

Finally, a trapezoidal rule is applied for the integration, i. e. 

riv-l 

3.3. Treatment of the boundary conditions. As we mentioned at the end of section 2, 
an appropriate discretization of boundary condition is important to exhibit very different 
phenomena. Therefore, we present respectively the numerical approximations for specular 
reflection boundary condition (2.6) and Neumann boundary condition (2.8). 

3.3.1. Numerical approximation for specular reflection boundary condition. The specular re- 
flection boundary condition in 2D reads as 

(3.6) 7^[/](^,x,v) = /(^,x,v'), 

with 

V • n = -v^ • n o v^ = v-2(vn)n 

where x g Fx = 50 is the point at the boundary, n is the interior normal at point x. We note 
that this specular reflection is just like a mirror. For example, if we follows the characteristic 
of the flux /(v), its reflected flow is corresponding to the velocity v^ 

We thus use a mirror procedure to construct / at each ghost point. For instance from the 
ghost point x^, we can find an inward normal n(xp), which crosses the boundary at x^ (see 
Figure 1). For velocity v, its reflected velocity with respect to x^ is = v-2(v-n(xp))n(xp). 
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Thus instead of computing / at the ghost point x^, we compute / at mirror point with respect 
to the boundary = 2xp x^ as fohows 

(3.7) /(x^,v) = /(x^,v'). 



iy + 2 
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Figure 1. Spatiahy two-dimensional Cartesian mesh. • is interior point, ■ is 
ghost point, □ is the point at the boundary, O is the point for extrapolation, 
the dashed line is the boundary. 



The last step is to approximate /(x^,v^) using / of interior domain. Let us assume that 
the values of the distribution function / on the grid points in Vt are given. We first construct 
a stencil E composed of grid points of for interpolation or extrapolation. For instance as it 
is shown in Figure 1, the inward normal n(xp) intersects the grid lines y = yiy^ Viy+i^ yiy+2 at 
points Pq, Pi , P2 . Then we choose the three nearest points of the cross point P^'^, / = 0, 1, 2, 
in each line, i.e. marked by a large circle. From these nine points, we can build a Lagrange 
polynomial q2{^) ^ Q2(IR^). Therefore we evaluate the polynomial q2{^) at x^, and obtain 
an approximation of / at the mirror point. We distinguish two cases of mirror points. In 
the case that mirror point x^ is surrounded by interior points, we interpolate / at mirror 
point x^; otherwise a WENO type extrapolation can be used to prevent spurious oscillations, 
which will be detailed below. 

Note that in some cases, we can not find a stencil of nine interior points. For instance 
when the interior domain has small acute angle sharp, the normal n can not have three cross 
points P^ - 0, 1, 2 in interior domain, or we can not have three nearest points of the cross 
point Pf , I = 0, 1, 2, in each line. In such a case, we alternatively use a first degree polynomial 
gi(x) with a four points stencil or even a zero degree polynomial g'o(x) with an one point 
stencil. We can similarly construct the four points stencil or the one point stencil as above. 
Two-dimensional WENO type extrapolation. A WENO type extrapolation ] was 
developed to prevent oscillations and maintain accuracy, which is an extension of WENO 
scheme [ ]. The key point of WENO type extrapolation is to define smoothness indicators. 
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which is designed to help us choose automaticahy between the high order accuracy and the 
low order but more robust extrapolation. Moreover a slightly modified version of the method 
was given in [ ] , such that the smoothness indicators are invariant with respect to the scaling 
of /. We now describe this method in 2D case. 

The substencils Sr, r = 0,1,2 for extrapolation are chosen around the inward normal n 
such that we can construct Lagrange polynomial of degree r. For instance in Figure 1, the 
three substencils are respectively 

{ So = {{xi^^Viy)}, 

S2 - {{xi^-i^yiy)^ (xi^^yiy), (xi^+i^yiy)^ (xi^-i^yiy+i), 

{xi^,yiy+i), (xi^+i,yiy+i), {xi^,yiy+2)^ 5 1/^^+2)5 (^i^+s? ^2^+2)}- 

Once the substencils Sr are chosen, we could easily construct the Lagrange polynomials in 

r r 
m=0/=0 

satisfying 

<?r(x) = /(x), XG 5^. 
Then the WENO extrapolation has the form 

2 

(3.8) /(x) = Wrqr(:><^), x g Sr, 

where are the nonlinear weights, which are chosen to be 



Wr = 



with 

a. 



v^2 



where £ = 10~^, do = + ^U^i di = \/Ax^ + ^y^, d2 = I - do - di. The coefficients are 
the smoothness indicators determined by 

and for r > 1, either we take = when /(x) = 0, Vx g Sr, or we choose 

Pr = ^ \. E r |i^|l-|-Hi?V(x))^dx, r = l,2, 

where cr is a multi-index and K = [xp- Ax/2, + Ax/2] x [yp-Ay/2,yp-\-Ay/2] and the point 
Xp is given by (xp,yp). 

3.3.2. Numerical approximation for Neumann boundary condition. In section 3.2, we have 
seen that discrete Laplace operator (3.5) is the only non-locally in space term in (2.7). In 
Figure 2, we illustrate an example of a discrete Laplace operator near the boundary. We note 
that Ni^^iy-i is not known, since x^ = (x^^,y^^_i) is out of domain. To approximate at 
ghost point x^, we have to use the boundary condition (2.8). 

In fact, if we denote n := {rix^riy), then the boundary condition (2.8) reads 

^xdxN + UydyN = 0. 
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Using a centered difference formula, it yields 

jV(x^)-jV(x,) jV(x^)-jV(x,) 
Tlx — + riy — = 0, 

~ ym ~ Vg 

where x^ := (x^,?/^), x^ := {xm^Vm) is the mirror point of the ghost point x^ with respect to 
the boundary. The previous equation is equivalent to 

(3.9) iV(x^) = iV(x3). 

Therefore instead of computing N at the ghost point x^, we extrapolate N at the mirror 
point x^ := {xm^Vm)- As shown in Figure 2, a nine points stencil 5^2 is found, thus we have 

8 

(3.10) N^^,^y-l = N{^m) = E^ziV(x,), X, G ^2, 

where Wi is weight of extrapolation calculated, for instance, by Lagrange polynomial. Finally 
by injecting (3.10) into (3.5), we complete the discretization. 

Remark 3.1. It is important to emphasize that the numerical scheme for the internal domain 
is disconnected to the numerical procedure for the treatment of the boundary. Therefore the 
present method can be applied to any other numerical scheme. The main point is to preserve 
the order of accuracy and then to eventually increase the stencil outside the domain. 




Figure 2. Illustration of discretization of Laplace operator near the 
boundary on two-dimensional Cartesian mesh. • is interior point, ■ is ghost 
point, □ is the point at the boundary, O is the point for extrapolation, the 
dashed line is the boundary. 
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Table 1 . Summary of the values used in the simulation 



Parameter 


Value 


run speed 


or -1 

= loiim ' s 


mean tumbling frequency 


1 o -1 
^0 = 35 


modulation or tumbling frequency or nutrient 


Xn = 60% 


modulation of tumbling frequency of chemoattractant 


XS - ^U/o 


stiffness of the response functions 


II5n = 1/Sn ^ 20s 


space scale 


X = 1mm 


time scale 


t = iOs 


doubling time 


T2 = ln2/r^2h 


degradation rate of the chemoattractant 


a = 5 • 10~^mol ■ 


production rate of the chemoattractant 


b = 4- 10^ cell-^s-^ 


consumption rate of the nutrient by the bacteria 


c = 2-W-'^cell-^s-^ 


diffusion coefficient of the nutrient molecules 


DN = 8-10-^cm^-s-l 


diffusion coefficient of the chemoattractant molecules 


Ds = 8- IQ-^cm^ ■ s-1 



4. Numerical examples 

In this section, we present a large variety of tests in 2dx and one dimensional in velocity 
space showing the validation of mathematical model and the effectiveness of our numerical 
schemes. In Test 1, we consider cell aggregation to study only chemotactic motility [ ]. 
In Test 2, we consider wave propagation formed by cells reorientation and the presence of 
nutrients in a disc. In Test 3, we study the interaction of two traveling wave in a, U shape. 
Then in Test 4, we consider again the traveling wave but in a more wide U shape, which is 
similar like the effectiveness of centrifugal force. Finally in Test 5 we add the effect of cell 
division for long time behavior using new parameters. 

4.1. Test 1 : cell aggregation. Motile bacteria are able to interact with their environment 
by accumulating in regions of high concentrations of certain chemicals called attractants and 
avoiding others called repellents. Cell division is not required to generate these multi-cellular 
structures because they form on a much shorter time scale than the cell-doubling time [ ]. 
The process leading to formation of multi-cellular clusters can be qualitatively understood 
as follows: fluctuations in the local cell density produce local gradients of attractants. Cells 
respond by moving up these concentration gradients thus amplifying the initial spatial non- 
uniformities in the cell distribution and forming multi-cellular clusters. 

We thus use the model of Othmer-Dunbar-Alt [22] to mimic the cell aggregation as follows 

5i/ + v-Vx/= f r[5](v,v')/(t,x,v')dv'- f r[5](v',v)/(t,x,v)dv', 

where the tumbling kernel T[S'](v^, v) describes the frequency of reorientation ^ v 

r[5](v', v) = V5 {dt log(5) + v' • Vx \og{S)) , 

and the chemical signal is secreted by the cells, following the reaction-diffusion equation (2.7). 

We use a square domain of size [-0.25, 0.25]^, which is uniformly divided by a mesh size of 
rixxriy = 120 X 120. The velocity space belongs to the unit circle 5^, and is uniformly divided 
into riy = 64 parts. All the parameters are chosen as in [^^8] and listed in Table 1. The initial 
chemoattractant Sq is equal to 0, and the initial distribution function fo is independent of 
the velocity v 

, , . 100m , , , ,o. 

/o(x,v) = exp(-100|xp), 

TT 

where m is the total mass. 
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The evolution of the ceh aggregation is presented in Figure 3. We observe that the initial 
density is a Gaussian distribution centered at (0,0). At time t = 0.2 1, the density forms a 
volcano profile, which disappears soon and is becoming into an exponential function at time 
t = It. Finally at time t = lOt, a steady exponential function is formed. The last observation 
is similar with the sharp boundary of clusters in [ ]. 




t = 



t = 





-0.2 -0. 



t = 0.2t 



t = 0.2t 




t = It 



t = It 



t = 10t 
(b) 





-0.2 -0.1 



t = 




t = 0.2t 




-0.2 -0. 



t = It 




Figure 3. Test 1 : Time evolution of the cell aggregation: (a) cell density 
in domain, (h) section plot of cell density, (c) section plot in log-scale of cell 
density. 



Furthermore, we observe in Figure 4 that with the total mass equal to m = tt/IOO, 7r/50 and 
7r/25 respectively, the size of clusters is almost the same. This phenomena coincides with the 
conclusion in [ ] that the steady-state size of clusters is almost independent of the number 
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of cells comprising them. In fact, following [ ] the steady density has a form 

p(x) ~ po exp(-A|x|), when |x| oo, 

where po is the maximum density, A depends on xs- Notice that the typical size of a cluster, 
defined as the mean radius given by 

^1^1^^ /x |x|p(x)rfx ^ 1 
/^p(x)dx A 

is independent of the total mass. This computation shows that the size of a cluster only 
depends on the model parameters and is entirely independent of the initial mass condition. 




Figure 4. Test 1 : The steady-state size of clusters comparison for different 
total masses at time t = 10 1: (a) section plot of cell density, (b) section plot 
in log-scale of cell density. 



4.2. Test 2 : wave propagation in a disc. Suspended Escherichia coli bacteria swim in 
convection-free geometries such as capillaries or micro-channels, collectively migrate towards 
nutrient-rich regions, in the form of propagation concentration waves [x]. In this test, we 
are interested in the model (2.1)-(2.8) to study concentration waves attracted by nutrient. 
Moreover, we consider a disc geometry to verify the numerical discretization of boundary 
conditions in section 3.3. 

The computational domain in space is a square domain of size [-3, 3]^, which is uniformly 
divided by a mesh size of Ux x riy - 80 x 80. The disc is inside of the square with radius 
equal to 3. It is clear that the boundary is not located on grids. Thus to achieve interior 
high order scheme, some artificial values on ghost points behind the boundary are needed. 
These values are given by using the numerical method presented in section 3.3. Moreover the 
velocity space belongs to the unit circle 5^, and is uniformly divided into = 64 parts. All 
the parameters are chosen as in [ ] and listed in Table 1. The initial distribution function 
/o is given by a Gaussian function 

/o(x,v) = poexp(|xp), 

where po is constant. The initial chemoattractant S^o is equal to and the initial nutrient A/'o 
is homogeneous equal to 1. 

The time evolution of concentration wave in a disc is illustrated in Figure 5. In the first 
row of Figure 5, we observe that the initial Gaussian density is extending and forming a 
propagating wave. When the circle wave arrives at the boundary, all cell are reflected by 
the boundary and attracted by nutrient remained in the disc center (see the second row of 
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Figure 5). In the third row of Figure 5, the circle wave contracts back to disc center, and 
finahy the cehs concentrates at the disc center. We notice that ceh diffusion appears when 
circle wave goes back to disc center. In fact, this diffusion is due to the stiffness of the 
response functions in the tumbling kernel [28]. 









3-2-10123 
X-axis 

t= 


-3-2-10123 
X-axis 

t = 4t 


-3-2-10123 
X-axis 

t = 8t 
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X-axis 
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X-axis 


-3-2-10123 
X-axis 

t = 19t 




■ o 11 




3-2-10123 
X-axis 

t = 23t 
Figure 5. Test ^ 


-3-2-10123 -3-2-10123 
X-axis X-axis 

t = 30f t = 33t 

I : Time evolution of the cell density at different time. 



4.3. Test 3 : interaction of traveling waves in a. U shape. In this test, we focus on 
the influence of the reorientation on the shape. The simulations are compared to a particular 
experiment by courtesy of Axel Buguin, Institut Curie (see the second row of Figure 9). We 
consider a channel of U shape, with initially homogeneous nutrient injected in the channel. 
Two clusters of bacteria are then imposed at two extremities of the channel. These two 
clusters move along the channel and finally meet at the channel center (the top of U shape). 
We note that these two clusters keep their bar shape till their meeting. 

To perform the simulation, we consider a channel of U shape with channel width equal to 
1 included in a rectangle computational domain [0, 8] x [0, 6], which is covered by an uniform 
mesh of size Uxxriy = 80x60. The velocity space belongs to the unit circle S^, and is uniformly 
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divided into riy = 64 parts. The numerical parameters used in the simulations are listed in 
Table 1. The initial distribution function fo is given by a constant as follows 

f 0.25, ifO<7/<l, 

/o(x,v) = 

[ 0, else. 

The initial chemoattractant Sq is equal to and the initial nutrient A^o is homogeneous equal 
to 1. 

The numerical simulations are presented in the first row of Figure 9. In Figure 9(a), we 
observe that two bar shape clusters form and are moving to the channel center. In Figure 9(b), 
two clusters go forward along the half-ring channel and keep well their bar shape. Finally in 
Figure 9(c), two clusters meet at the channel center. We see that the numerical simulations 
have a good agreement with the experiment. 




(a) (b) (c) 

Figure 6. Test 3 : Time evolution of the cell density (a) t= 0.65 sec. (h) t 
= 2.65 sec. (c) t = 4-53 sec. At the top numerical results and at the bottom 
experiments on Escherichia coli (courtesy of Axel Buguin, Institut Curie). 



4.4. Test 4 : one traveling wave in a ?7 shape. In the last test, we consider again 
traveling wave in a ?7 shape but with a more wide channel than the previous one. The 
experience is shown in the second row of Figure 7 and Figure 8. Again we inject homogeneous 
nutrient in the channel, but we consider only one cluster of bacteria at the right extremity 
of the channel. We observe that at straight part of channel the cluster goes ahead in a bar 
shape. Once the cluster enters into the half-ring part, the bacteria near interior circle goes 
faster than the one near the exterior circle. Moreover, bacteria contracts towards the exterior 
circle. Before the cluster enters into the left straight part of channel, bacteria concentrate 
almost near the exterior circle. When the cluster goes forward the other extremity, the cluster 
recovers its original bar shape. 

To perform the simulation, we consider a channel of U shape with channel width equal 
to 3 included in a rectangle computational domain [0,6.5] x [0,8], which is covered by an 
uniform mesh of size n^; x = 65 x 80. The velocity space belongs to the unit circle 5^, and 
is uniformly divided into = 64 parts. The numerical parameters used in the simulations 
are listed in Table 1. The initial condition is the same as in the previous test. 
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The numerical simulations are presented in the first row of Figure 7 and Figure 8. We 
can see that the numerical simulations are very similar as the experience one. In fact, this 
phenomena is due to the directional persistence of chemotactic bacteria in a traveling con- 
centration wave. When bacteria enter into a wide half-ring, they keep going straight and 
accumulate by the reflection of the exterior circle. It is very similar like the effectiveness of 
centrifugal force, and can be observed significantly in a wide channel. This test shows that 
the model (2.1)- (2.8) represents well the chemotactic bacteria behavior and our numerical 
discretization based on Cartesian mesh approximates well the continuous model. 
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Figure 7. Test 4 : Numerical simulations (top) and experiments on Es- 
cherichia coli (bottom) : time evolution of the cell density (a) t = 7.5 t (b) t 
= 18.5 t (c) t = 23tin sec. 



4.5. Test 5 : long time behavior and pattern formations. In this last numerical test, 
we consider the full kinetic model with cell divisions and degradation 



(4.1) 



dl 
dt 



+ v-Vx/ = Q(/) + r/-7e(p)e(p-poo), xef2, ^reV, 



where the division rate is given by [26] 



GoN 



it takes into account two experimental facts : the slowing down of the growth rate for low 
nutrients concentrations and its finite quantity for high concentrations. The third term on the 
right hand side of (4.1) describes the vegetative cells into anabiotic form due to the increase 
of the local density. The transition starts when the total cell density reaches the value poo 
and O is the Heaviside function. 




(d) (e) (f) 

Figure 8. Test 4 : Numerical simulations (top) and experiments on Es- 
cherichia coli (bottom) : time evolution of the cell density (d) t = 28 t (e) t 
= 34.5 i (f) t = 46.5 t in sec. 



Actually such a source term has been proposed in [ ] in the framework of a macroscopic 
Patlak-Keller-Segel model where the unknown is the total density p. Numerical simulations 
of this macroscopic model have shown pattern formations as the ones observed in experiments 
[26, 6]. Here we consider the following initial density 



/o(x,v) := I J 



1 if|x|<l, 
else 



and the initial nutrient concentration is uniform A^o = 0.5. In the nutrient and chemoattrac- 
tant system of equations (2.7) we take the following parameters Ds = = 1, the production 
rate of chemoattractant 6 = 20, the degradation rate of chemoattractant a = 8 and the con- 
sumption rate of nutrient c = 0.8, whereas in (4.1), we choose a = 0.1, poo = 15 and for the 
turning operator (2.2)-(2.4),we have S = 20, i/jq = 1, xa^ = 1/2 and xs = 1/10. The numerical 
simulations are presented in Fig. 9. We observe that for high initial nutrient concentration, 
cell density in the expanding ring becomes sufficient both for the break of stability of the 
uniform cell distribution and for their aggregation. In particular, if after the formation of the 
successive set of aggregates the expanding ring had time to grasp a certain part of the cells, 
participating in aggregation, then a radial pattern is formed. 



5. Conclusion 

In this paper we present a new algorithm based on a Cartesian mesh for the numerical ap- 
proximation of kinetic models on an arbitrary geometry boundary modelling chemosensitive 
movements. We present first a kinetic model for chemotactic bacteria interacting with two 
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Figure 9. Test 5 : Time evolution of the cell density (a) t = t (h) t = 15 
t (c) t = 25t (d) t = 35t(e)t = 45t and (f) t = 55tin sec. 



chemical substances, i.e. nutrient and chemottractant. Then we give the numerical discretiza- 
tion for this kinetic model and the numerical method for the boundary conditions based on a 
Cartesian mesh. A large various numerical tests in 2Dx x are shown and compared with 
biological experiences. We conclude that on the one hand this kinetic model represents well 
the chemotactic bacteria behaviors and on the other hand our numerical method is accurate 
and efficient for numerical simulation. 
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